function [IRFsdf,IRFsdfnoEZ,IRFsdfnoD] = computeIRFsdf(model,shockSize,xf,xs,appOrder,IRFlength)

invI_hx  = (eye(model.nx)-model.hx)\eye(model.nx);
if isempty(xf)
    xf = invI_hx*model.h0;
end

rng(1,'twister');
numSim = 1000;
shocksMat                     = zeros(size(model.eta,2),IRFlength,numSim);
shocksMat(:,:,1:numSim/2)     = randn(size(model.eta,2),IRFlength,numSim/2);
shocksMat(:,:,numSim/2+1:end) = -shocksMat(:,:,1:numSim/2);

% We compute several sample paths with the shocks and without the shock
nx    = model.nx;
nx1   = model.nx1;
ne    = model.ne;
heteroShocks = reshape(model.heteroShocks,model.ne,1);
xSim_pshock  = zeros(model.nx,IRFlength,numSim);
ySim_pshock  = zeros(model.ny,IRFlength,numSim);
xSim         = zeros(model.nx,IRFlength,numSim);
ySim         = zeros(model.ny,IRFlength,numSim);
for s = 1:numSim
    %% WIth a shock
    xf_cu = xf;
    xs_cu = xs;
    for t=1:IRFlength
        xf_cup = model.h0 + model.hx*xf_cu;       
        if t == 1
            xf_cup(nx1+1:nx,1) = xf_cup(nx1+1:nx,1) +...
                (2./(exp(-heteroShocks(:).*xf_cu(nx1+1:nx,1))+1)).*model.eta(nx1+1:nx,:)*(shockSize+shocksMat(:,t,s));
        else
            xf_cup(nx1+1:nx,1) = xf_cup(nx1+1:nx,1) +...
                (2./(exp(-heteroShocks(:).*xf_cu(nx1+1:nx,1))+1)).*model.eta(nx1+1:nx,:)*shocksMat(:,t,s);
        end
        if model.appOrder > 1
            AA_cu  = kron3(xf_cu,xf_cu);
            xs_cup = model.hx*xs_cu + model.Hxx*AA_cu;
        end
        
        if model.appOrder == 1
            ySim_pshock(:,t,s) = model.g0 + model.gx*xf_cup;
            xSim_pshock(:,t,s) = xf_cup;
        elseif model.appOrder == 2
            AA_cup      = kron3(xf_cup,xf_cup);
            ySim_pshock(:,t,s) = model.g0 + model.gx*(xf_cup+xs_cup) + model.Gxx*AA_cup;
            xSim_pshock(:,t,s) = xf_cup+xs_cup;
        end
        % Updating the states (dev from steay state)
        xf_cu = xf_cup;
        if model.appOrder > 1
            xs_cu  = xs_cup;
        end
    end

    %% The benchmark without a shock
    xf_cu = xf;
    xs_cu = xs;
    for t=1:IRFlength
        xf_cup = model.h0 + model.hx*xf_cu;       
        xf_cup(nx1+1:nx,1) = xf_cup(nx1+1:nx,1) +...
               (2./(exp(-heteroShocks(:).*xf_cu(nx1+1:nx,1))+1)).*model.eta(nx1+1:nx,:)*shocksMat(:,t,s);
        if model.appOrder > 1
            AA_cu  = kron3(xf_cu,xf_cu);
            xs_cup = model.hx*xs_cu + model.Hxx*AA_cu;
        end        
        if model.appOrder == 1
            ySim(:,t,s) = model.g0 + model.gx*xf_cup;
            xSim(:,t,s) = xf_cup;
        elseif model.appOrder == 2
            AA_cup      = kron3(xf_cup,xf_cup);
            ySim(:,t,s) = model.g0 + model.gx*(xf_cup+xs_cup) + model.Gxx*AA_cup;
            xSim(:,t,s) = xf_cup+xs_cup;
        end
        % Updating the states (dev from steay state)
        xf_cu = xf_cup;
        if model.appOrder > 1
            xs_cu  = xs_cup;
        end
    end
end
%% For the SDF
% Recall 
%mReal_cup/pai_cup = BETTA*d_cup/d_cu*...
%    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
%    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0)/pai_cup;  
% so 
%mReal_cu = BETTA*d_cu/d_ba1*...
%    (AA*evf_ba1^(1/(1-ALFA))/(vf_cu*muz_cu^((1-CHI)*(1-CHI0)) ))^ALFA...   
%    *(c_cu-B*c_ba1/muz_cu)^(-CHI)/((c_ba1-B*c_ba2/muz_ba1)^(-CHI))*muz_cu^(-CHI*(1-CHI0)-CHI0)/pai_cu;  

% Position of variables
posC      = find(strcmp(model.labely,'$c_t$'));
posVF     = find(strcmp(model.labely,'$vf_t$'));
posEVF    = find(strcmp(model.labely,'$evf_t$'));
posPAI    = find(strcmp(model.labely,'$\pi_t$'));

posD      = find(strcmp(model.labelx,'$d_{t}$'));
posmuz    = find(strcmp(model.labelx,'$\mu _{z,t}$'));
posC_ba1  = find(strcmp(model.labelx,'$c_{t-1}$'));

% Get model parameters
BETTA = model.params.BETTA;
B     = model.params.B;
CHI   = model.params.CHI;
CHI0  = model.params.CHI0;
ALFA  = model.params.ALFA;
AA    = model.params.AA;

% Value of the states and controls at time 0
if model.appOrder == 1
    y0 = model.g0 + model.gx*xf;
    x0 = xf;
elseif model.appOrder == 2
    y0    = model.g0 + model.gx*(xf+xs) + model.Gxx*kron3(xf,xf);
    x0    = xf_cu+xs_cu;
end

sdf_pshock = zeros(IRFlength,numSim);
sdf        = zeros(IRFlength,numSim);
sdfnoEZ_pshock = zeros(IRFlength,numSim);
sdfnoEZ        = zeros(IRFlength,numSim);
sdfnoD_pshock  = zeros(IRFlength,numSim);
sdfnoD         = zeros(IRFlength,numSim);
for s = 1:numSim
    %% WIth a shock
    for t=1:IRFlength
        % SDF_cu = BETTA*d_cu/d_ba1*...
        %    (AA*evf_ba1^(1/(1-ALFA))/(vf_cu*muz_cu^((1-CHI)*(1-CHI0)) ))^ALFA...   
        %    *(c_cu-B*c_ba1/muz_cu)^(-CHI)/((c_ba1-B*c_ba2/muz_cu)^(-CHI))*muz_cu^(-CHI*(1-CHI0)-CHI0)/pai_cu;  
        if t == 1
            c_ba2   = x0(posC_ba1,1);
            c_ba1   = xSim_pshock(posC_ba1,t,s);
            d_cu    = xSim_pshock(posD,t,s);
            d_ba1   = x0(posD,1);
            muz_cu  = xSim_pshock(posmuz,t,s);
            muz_ba1 = x0(posmuz,1);
            
            c_cu    = ySim_pshock(posC,t,s);
            evf_ba1 = y0(posEVF,1);
            vf_cu   = ySim_pshock(posVF,t,s);
            pai_cu  = ySim_pshock(posPAI,t,s);
        else
            c_ba2   = xSim_pshock(posC_ba1,t-1,s);
            c_ba1   = xSim_pshock(posC_ba1,t,s);
            d_cu    = xSim_pshock(posD,t,s);
            d_ba1   = xSim_pshock(posD,t-1,s);
            muz_cu  = xSim_pshock(posmuz,t,s);
            muz_ba1 = xSim_pshock(posmuz,t-1,s);

            c_cu    = ySim_pshock(posC,t,s);
            evf_ba1 = ySim_pshock(posEVF,t-1,s);
            vf_cu   = ySim_pshock(posVF,t,s);
            pai_cu  = ySim_pshock(posPAI,t,s);
        end
        sdf_pshock(t,s) = BETTA*exp(d_cu)/exp(d_ba1)*...
           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           *(exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  
        
        sdfnoEZ_pshock(t,s) = BETTA*exp(d_cu)/exp(d_ba1)*...           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           (exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  

        sdfnoD_pshock(t,s) = BETTA*...   *exp(d_cu)/exp(d_ba1)*...
           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           *(exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  
    end


    %% The benchmark without a shock
    for t=1:IRFlength
        if t == 1
            c_ba2   = x0(posC_ba1,1);
            c_ba1   = xSim(posC_ba1,t,s);
            d_cu    = xSim(posD,t,s);
            d_ba1   = x0(posD,1);
            muz_cu  = xSim(posmuz,t,s);
            muz_ba1 = x0(posmuz,1);
            
            c_cu    = ySim(posC,t,s);
            evf_ba1 = y0(posEVF,1);
            vf_cu   = ySim(posVF,t,s);
            pai_cu  = ySim(posPAI,t,s);
        else
            c_ba2   = xSim(posC_ba1,t-1,s);
            c_ba1   = xSim(posC_ba1,t,s);
            d_cu    = xSim(posD,t,s);
            d_ba1   = xSim(posD,t-1,s);
            muz_cu  = xSim(posmuz,t,s);
            muz_ba1 = xSim(posmuz,t-1,s);

            c_cu    = ySim(posC,t,s);
            evf_ba1 = ySim(posEVF,t-1,s);
            vf_cu   = ySim(posVF,t,s);
            pai_cu  = ySim(posPAI,t,s);
        end
        sdf(t,s) = BETTA*exp(d_cu)/exp(d_ba1)*...
           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           *(exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  

        sdfnoEZ(t,s) = BETTA*exp(d_cu)/exp(d_ba1)*...           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           (exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  

        sdfnoD(t,s) = BETTA*...   *exp(d_cu)/exp(d_ba1)*...
           (AA*exp(evf_ba1)^(1/(1-ALFA))/(exp(vf_cu)*exp(muz_cu)^((1-CHI)*(1-CHI0)) ))^ALFA...   
           *(exp(c_cu)-B*exp(c_ba1)/exp(muz_cu))^(-CHI)/((exp(c_ba1)-B*exp(c_ba2)...
           /exp(muz_ba1))^(-CHI))*exp(muz_cu)^(-CHI*(1-CHI0)-CHI0)/exp(pai_cu);  
    end
end
IRFsdf     = mean(sdf_pshock,2)-mean(sdf,2);
IRFsdfnoEZ = mean(sdfnoEZ_pshock,2)-mean(sdfnoEZ,2);
IRFsdfnoD  = mean(sdfnoD_pshock,2)-mean(sdfnoD,2);
